/ 



APPLICATION 
FOR 

UNITED STATES LETTERS PATENT 



APPLICANT NAME: Ronald L Gordon 
TITLE: METHOD FOR COMPUTING PARTIALLY COHERENT AERIAL IMAGERY 

DOCKET NO.: FIS9-2003-0254US1 

INTERNATIONAL BUSINESS MACHINES CORPORATION 
NEW ORCHARD ROAD, ARMONK, NY 10504 



CERTIFICATE OF MAILING UNDER 37 CFR 1.10 
I HEREBY CERTIFY THAT, ON THE DATE SHOWN BELOW 
THIS CORRESPONDENCE IS BEING DEPOSITED WITH 
THE UNITED STATES POSTAL SERVICE IN AN 
ENVELOPE ADDRESSED TO THE U.S. PATENT AND 
TRADEMARK OFFICE, ARLINGTON, VA 22202 AS 
"EXPRESS MAIL POST OFFICE ADDRESSEE- 
MAILING LABEL # EU702258644US 



On 



07-17-03 



Karen Cinq-Mars 
Name of person 




Signature 



ig paper 



Date 



METHOD FOR COMPUTING PARTIALLY COHERENT AERIAL IMAGERY 
TECHNICAL FIELD 

The present invention relates in general to manufacturing processes that require 
lithography and, in particular, to methods of designing photomasks and optimizing 
5 lithographic and etch processes used in microelectronics manufacturing. 

BACKGROUND INFORMATION 

During microelectronics manufacturing, a semiconductor wafer is processed 
through a series of tools that perform lithographic processing, followed by etch 

10 processing, to form features and devices in the substrate of the wafer. Such processing 
has a broad range of industrial applications, including the manufacture of 
semiconductors, flat-panel displays, micromachines, and disk heads. 

The lithographic process allows for a mask or reticle pattern to be transferred via 
spatially modulated light (the aerial image) to a photoresist (hereinafter, also referred to 

1 5 interchangeably as resist) film on a substrate. Those segments of the absorbed aerial 

image, whose energy (so-called actinic energy) exceeds a threshold energy of chemical 
bonds in the photoactive component (PAC) of the photoresist material, create a latent 
image in the resist. In some resist systems the latent image is formed directly by the 
PAC; in others (so-called acid catalyzed photoresists), the photo-chemical interaction 

20 first generates acids which react with other photoresist components during a post- 

exposure bake to form the latent image. In either case, the latent image marks the volume 
of resist material that either is removed during the development process (in the case of 
positive photoresist) or remains after development (in the case of negative photoresist) to 
create a three-dimensional pattern in the resist film. In subsequent etch processing, the 

25 resulting resist film pattern is used to transfer the patterned openings in the resist to form 

an etched pattern in the underlying substrate. 

Diffraction, interference and processing effects that occur during the transfer of 
the image pattern causes the image or pattern formed at the substrate to deviate from the 
desired (i.e. designed) dimensions and shapes. These deviations depend on the 
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interaction of the pattern configurations with the process conditions, and can affect the 
yield and performance of the resulting microelectronic devices. Various techniques have 
been used to compensate for and correct for these deviations. Such techniques include 
known as optical proximity correction (OPC), for example, by biasing selected mask 

5 features to compensate for deviations. Other techniques include using sub-resolution 

assist features (SRAFs), also known as scattering bars or intensity leveling bars, to 
improve the uniformity of grating characteristics of the mask, and thereby improve 
optimization of lithographic process conditions for the mask. Phase shifted mask 
technology (PSM) has also been used to improve the contrast of image features by 

1 0 destructive interference, and thus improve resolution. These and other various techniques 
for improving the lithographic process are generally referred to as resolution 
enhancement techniques (RETs). 

Prediction of the partially coherent images resulting from illumination in a 
modern lithographic scanner is of paramount importance as looming technology nodes 

15 stress the use of Resolution Enhancement Technologies (RETs), such as SubResolution 

Assist Features (SRAFs) and Alternating Phase-Shift Masks (AltPSM). Without fast and 
accurate simulation, it would be impossible to employ an strong RET solution in a 
practical setting. The reason for this is because simulations have made the transition 
from learning/research tool to a major ingredient in the design stage. However, in 

20 practice, only a portion of a mask pattern can be simulated at a time, to allow for 
reasonable computation times. 

The ability to accurately predict the resulting aerial image, latent image and/or 
etched pattern due to the lithographic and etch processes is crucial for ensuring sufficient 
manufacturing yields and reducing costs of manufacturing. The aerial image of a mask 

25 pattern is the distribution of intensity at the plane of the wafer or resist surface. The 

accurate simulation of the aerial image is key in the design of photomasks, for example, 
by model-based optical proximity correction (model-based OPC). In model-based OPC, 
for given lithographic process conditions (e.g. illumination source parameters, numerical 
aperture (NA) and partial coherence (a 0 ) of the lithography tool, specified exposure dose 
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and defocus, etc.), and an initial mask design, the resulting image is simulated. In the 
model-based OPC process, the simulated image is compared to the target (desired) 
design, and deviations from the desired image are determined. The mask design is 
iteratively modified until the simulated image agrees with the target design within an 
5 acceptable tolerance. Thus the accuracy of the simulated image is crucial in obtaining 

viable mask designs, and the speed of such calculations impacts the cost of designing the 
masks. 

Simulated images can also be used to improve understanding the interaction of, 
and for optimizing the lithographic process conditions to maximize resolution and 

10 improve yields. For example, for a particular design pattern, a decision to use one type of 
resolution enhancement technique over another, for example, whether to use altPSM or 
SRAFs involves an understanding what the best range of process windows will be for a 
range of altPSM or SRAF processes. A wide range of process conditions must be 
simulated and corresponding images simulated for each process condition, at a required 

15 accuracy. The simulation of a single aerial image using conventional methods, with 

sufficient accuracy (e.g. within 1%), typically takes hours or even days. 

To get an idea of what sort of accuracy is required in simulation, consider that the 
aerial image simulator can be viewed as a metrology tool. There is an inherent 
uncertainty in metrology, for example in using metrology techniques such as (scanning 

20 electron microscope) SEM, in which a target to be measured is bombarded with 

electrons, which in turn produce a signal, indicating, for example, line widths in the 
target. However, there is an inherent uncertainty in these measurements, for example, 
caused by charge damage to features on the target. The precision to tolerance ratio (P to 
T ratio) is the ratio of the precision, or accuracy, of the metrology tool to the tolerance 

25 specification for the device being measured. The specification for line widths (CD) may 

be, for example, 90 nm, within 3 sigma. The demand on the line width distribution is 
such that CD has a mean value of 90 nm, wherein the 3 sigma variation is within 5% of 
90 nm (e.g. ±4.5 nm tolerance). A typical spec for P-T ratio is to measure the line within 
a small fraction of 4.5 nm (e.g. 20% or 0.90 nm or 9.0 Angstroms). A simulation tool 
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should provide numerical accuracy in a similar vein as the metrology tools specifications, 
e.g. within .90 nm. One conventional method of simulating aerial images is to use a 
gridding algorithm, as in the prior art outlined above, but in order to obtain the precisions 
required, to obtain the required precision, the smaller grid sizes result in a large number 
5 of gridding intervals, which in turn result in impractical computation times. Such 
gridding methods cannot be used to simulate large portions of a mask in a practical 
amount of time. 

Partially coherent imagery is simulated using what is now called the Hopkins 
Model (see, H.H. Hopkins, "On the diffraction theory of optical images," Proc. of the 

10 Royal Society of London, Vol. A 217, pp. 408-432 (1953)). The Hopkins model for 
simulating the aerial image is a method to compute aerial images using the frequency 
space information of the optical system. In doing this, the calculation can be split into 
two independent steps. Here, the intensity, defined as the average of the square 
magnitudes of the electric fields that emanate from each point in the illumination source, 

1 5 is expressed as a quadratic form in the mask spectrum - which is the Fourier transform 
(in the spatial frequency domain) of the mask transmittance - multiplied by a matrix of 
complex numbers called Transmission Cross-Coefficients (TCCs), which takes into 
account the source (illumination shape) and pupil (aberrations, defocus, vector, obliquity) 
information. The TCCs are autocorrelations of the transfer function of the pupil 

20 (henceforth called the "pupil function") weighted by the spatial Fourier transform of the 
illumination source (referred to hereinafter as the "source function"). The pupil is the 
image of the limiting stop, or aperture, in an optical system, but for the present purposes, 
the entrance/exit pupils represent the input/output planes for the optical system. These 
autocorrelations are, in general, double integrals over a very complicated region defined 

25 by the intersection of the pair of frequency-offset pupils and the source. Because of the 
complexity of the shape of this region, computation of these TCCs is potentially 
expensive. 

The second step is the computation of the mask spectrum (or, Fourier transform of 
the mask, referred to hereinafter as the "mask function"). This can be computed 
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analytically for most, if not all, one-dimensional masks in lithography. An analytical 
computation of the Fourier transform for all two-dimensional shapes can be obtained, for 
example, by the method disclosed in U.S. Patent Application Ser. No. 10/353,900, which 
is commonly assigned to the Assignee of the present application. This analytical 
5 calculation, while exact, tends to be expensive compared with the decomposition 

techniques in the SOCS method, because each edge required a trigonometric function 
evaluation. Nevertheless, such a calculation leaves little doubt as to the overall accuracy 
of the computation. 

A major simplification is typically made by assuming that the portion of the mask 
10 to be simulated is periodic. This is because the spectrum of periodic gratings is nonzero 
at a discrete number of frequencies, meaning that the frequencies that have nonzero TCCs 
lie on a regular, discrete grid in spatial frequency space. This periodic assumption is used 
by all current lithography simulators, since rigorous image computation over an entire 
(26 mm) field is impractical. This, of course, is an approximation in itself, as real 
1 5 circuitry to be manufactured will not necessarily possess this periodic property. 

Nevertheless, a single local set of shapes in a region of interest, or cell, of any size can be 
simulated with great accuracy from the periodic assumption. The reason for this is that 
optical diffraction phenomena have a range of CZ/(NAcr 0 ), where NA is the numerical 

aperture, X is the wavelength, cr 0 is the partial coherence (or pupil-fill), and C is a 

20 constant equal to about at most 10 (for example, see Wong et al., U.S. Patent 6,223,139). 
Therefore, there is a finite range of influence, beyond which features have no impact on 
the desired pattern. Thus, if ones assumes periodicity (or any other boundary condition) 
for a unit cell of the nominal pattern size plus this influence range, then the simulated 
image inside the unit cell can be considered exact for all practical purposes. 

25 With the mask periodic and the spectrum discrete, the set of TCCs to be evaluated 

is therefore discrete and finite in number. In general, they fill a four-dimensional matrix, 
and it can be shown that the number of TCCs needed for computation varies as the square 
of the area of the unit cell - in other words, as the fourth power of the length of a nearly 
square cell. Therefore, even with this simplification, the number of difficult, double 
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integrals that are needed to accurately define the image becomes unmanageable for even 
moderately large (< 10 jam 2 area) cells. 

In all current lithography simulation software, the TCCs are currently computed 
by gridding the source and pupil for a numerical integration. This gridding can be 

5 sophisticated: for example, in the software SPLAT (see K.K.H. Toh and A. Neureuther, 

"Identifying and Monitoring Effects of Lens Aberrations in Projection Printing," SPIE, 
Vol. 772, pp. 202-209, Optical Microlithography VI: 4-5 March 1987, Santa Clara, CA), 
a fixed grid is used in one direction, while an adaptive grid (it refines until a certain 
tolerance is reached) is used in another. In order to achieve this, the integration limits of 

10 the integration region must be computed. This allows for greater accuracy in computing 
aberrated images, for example. Unfortunately, because of the adaptive stepping in the 
integration, the algorithm runs as long as it needs in order to achieve a certain accuracy. 
This can take a long time, especially with pupils that have large phase variations, such as 
in large defocus and/or aberrations. 

1 5 The problem of providing fast simulation has already been attacked though the 

use of the so-called "Sum of Coherent Systems" (SOCS) method detailed in N. Cobb and 
A. Zakhor, "Fast low-complexity mask design", Proc. SPIE, Vol. 3334, pp. 313-327 
(1995). This method expresses the intensity as a sum over the squares of convolutions of 
the mask transmittance function with coherent point spread functions, or kernels; these 

20 kernels are inverse Fourier transforms of the eigenfunctions of the TCC matrix. By 
expressing the mask in terms of a discrete set of elements, these convolutions can be 
precomputed and stored in a database. This allows for very rapid image computation that 
is crucial for model-based OPC (MBOPC). There are caveats with this SOCS 
methodology, however. First of all, the sum used must be finite (and is typically between 

25 8 and 12 terms), whereas the exact expression involves an infinite number of terms; 

typical accuracy estimates from using the finite series for the aerial image intensity are 
between 80 and 90%, with a relatively slow rate of convergence. While this may suffice 
for past generations with a larger error tolerance (that is, the larger allowed linewidth 
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variation could accept larger simulation errors, as noted above), this will not serve the 
future generations, where the error tolerances are rapidly vanishing. 

Second, these kernels, are eigenfunctions of a matrix operator defined by the 
TCCs. In order to compute these eigenfunctions, the TCCs must be computed. Since the 
5 TCCs are independent of the mask function, these TCCs will be computed only once, 
whereas the intensity is computed thousands of times. This is true if only the mask 
features (and not the parameters of the projection system) are to be varied during a 
MBOPC session. It turns out, however, that the task of calibration of a model to data has 
gotten so difficult that some physical parameters that are known to be in error - such as 
10 the partial coherence or focus - are varied as well. The fitting of these physical 

parameters demands recomputation of the TCCs, so the ability to compute them very 
quickly is crucial here. 

There is a method of computing the TCC matrix rapidly due to Liebchen (U.S. 
Patent Application Pub. No. US 2002/0062206 Al). Here, the TCC matrix itself is 

1 5 approximated as a bilinear form rj Atj , where rj is a vector of orthonormal Hermite- 

Gaussian terms, and A is a Hermitian matrix. The benefit of doing this is that the size of 
the TCC matrix can be reduced substantially, making changes with respect to optical 
parameters such as the numerical aperture or defocus very fast. While this certainly has 
resulted in a significant increase in computational speed, it is not obvious what accuracy 

20 of the image intensity that this method produces, as not only the pupil and source are 

evaluated on a grid and numerical integration is necessary on at least one dimension, but 
the expression of the TCC matrix in this orthogonal basis is only exact when an infinite 
number of terms are used. 

Accuracy is crucial, and although few have ever questioned the half-century-old 

25 theory of optical image formation via partially coherent Kohler illumination (despite the 

scant experimental verifications of the various parametric dependencies, especially those 
relating to "obliquity factors" which are generally considered to have small impact), 
numerical accuracy in proprietary software has always been a concern. Many 
benchmarking projects have been performed, but most have been simple numerical 
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checks between software. Speed is also a major issue, because many thousands of 
simulations must take place in order to perform a single optimization. Historically, there 
has always been a tradeoff between speed and accuracy, in that speed requires one to 
sample less source integration points, or less eigenfunction kernels in the "fast" 
5 algorithms used in Model-Based OPC software. On the other hand, analytical solutions, 

if they exist, solve the speed/accuracy issue, since all the computation time comes from 
the numerous double integrations required for the image computation. If one were to 
replace the double integrations (thousands of function evaluations) with a single function 
evaluation, then the speed of the algorithm would increase many times over, and the 
10 accuracy would be to machine precision. The trouble is, of course, finding such a 
solution. 

Rigorous analytical methods of aerial images would be preferred to the numerical 
gridding methods commonly in use today, but are typically impractical because the 
resulting integrals cannot directly be evaluated in many cases. However, such analytical 

15 methods have been available for computing aerial images in special cases. The simplest 
case of this was illustrated by Kintner (Kintner, Eric C, "Method for the calculation of 
partially coherent imagery," Applied Optics, Vol. 17, No. 17, pp. 2747-2757, Sept. 1, 
1978) for simulation of in-focus, partially coherent images of one-dimensional (ID) 
gratings. By considering all of the possible geometrical configurations of the TCCs in 

20 ID (i.e., along the x-axis), Kintner was able to compute the exact value of the TCCs due 

to the fact that these values are equal to the area of a region bounded by 3 circles, all of 
whose centers are collinear. While this was a useful first step, the important cases of 
more general patterns seen in lithography were untreated, as well as the ability to 
compute defocus effects. 

25 Going one step further, Subramanian (Subramanian, S., "Rapid calculation of 

defocused partially coherent images," Applied Optics, Vol. 20, No. 10, pp. 1854-1857, 
May 15, 1981) considered the simulation of defocused images of the one-dimensional 
(ID) gratings. Here, Kintner' s work was built upon by finding the integration limits 
within each geometrical configuration. In the paraxial (i.e., small numerical aperture, 
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with NA < 0.4) approximation, the defocus term could be integrated analytically in one 
dimension, while numerical integration was used for the other. While this extended the 
art significantly, it is still limited to one-dimensional (ID) gratings and small values of 
the numerical aperture. 

The method outlined by Liebchen takes the approach of decomposing these TCCs 
(i.e. representing the TCC matrix as a bilinear Hermite-Gaussian expansion). While 
Liebchen's method does potentially reduce the number of TCCs needed, it also 
introduces a gridding; the combination of the decomposition and gridding methodologies 
introduce potential inaccuracies. Further, it is unclear how large values of defocus or 
aberrations are treated with any accuracy, as all treatments are Taylor expansions of the 
phases in the frequency variable. For small aberrations, this is acceptable, but such 
expansions quickly lose accuracy in the face of even a moderate defocus 

Therefore, there is a need to compute the TCCs without recourse to grids or 
orthogonal basis decompositions so that accuracy is maintained, or at least monitored, 
while the speed needed to computed moderately large patterns is satisfied. 
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SUMMARY OF THE INVENTION 

It is an objective of the present invention to provide a method for obtaining solely 
analytical expressions for all of the TCC integrals. A further objective of the present 
invention is to provide an approximate analytical representation for the TCCs that is 

5 accurate to within a desired small error, for example, within the precision of a given 
computer, which is typically on the order of about 10" 6 for typical single precision 
calculations in current machines. It is yet a further objective of the present invention to 
provide a method for computing aerial images very rapidly compared to conventional 
numerical integration techniques, and to within a desired accuracy, for example very 

10 close to machine precision. That is, a primary objective of the present invention is to 
eliminate the tradeoff between speed and accuracy in the simulation of aerial images, 
achieving both simultaneously. 

In accordance with the present invention, an analytical solution for the TCC 
integral has been derived. For the in-focus, aberration-free case, the analytical solution is 

1 5 exact, using only a single arctangent and at most 4 square roots per transmission cross- 

coefficient (TCC). This solution has been extended to defocus and small aberrations, and 
has been implemented in computer code. 

The invention takes advantage of two facts concerning the TCCs: 1) The 
integration region of the TCCs is bounded by circular arcs for circular-shaped sources 

20 and pupils; and 2) the integrand of the TCC expression can be represented by a simple, 
linear phase for the case of a defocused pupil in the paraxial approximation. These facts, 
combined with the application of Stokes' Theorem (a basic theorem of integral calculus) 
results in a reduction of the double integral to a single integral for any arbitrary 2D 
pattern. 

25 In a further aspect of the present invention, the TCC integration regions are 

characterized all in terms of a finite set of all the possible geometrical configurations, 
based on a rotational alignment of the spatial coordinate axis with the axis of symmetry 
of the integration region. These regions change shape as the spatial frequencies vary, as 
do the integration limits over which the integrand is to be evaluated. It turns out that 
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there are 18 distinct geometrical configurations when the double integral is used. 
Furthermore, when Stokes' Theorem is applied, only the boundary of these geometric 
configurations are considered, and it turns out that the number of distinct configurations 
is reduced to 9. Such a reduction simplifies the program logic and speeds up 
5 calculations. 

In a further aspect of the present invention, all of the integrals over the various 
arcs that comprise the boundaries of the finite geometrical regions is reduced to the same 
integral form. Since this particular integral does not have an exact representation in 
terms of purely analytical functions, the present invention provides an approximation for 

10 the integrand that differs from the original function in absolute value by a desired small 
number, for example, on the order of single precision over all possible values of the 
integration limits. In accordance with the present invention, this new function is chosen 
such that the error substantially monotonically decreases as the number of terms used in 
the approximation increases. Furthermore, the integral of the terms of the approximate 

15 function is representable in terms of analytical functions that can be computed with a few 

arithmetic steps. Therefore, the present invention permits computation of a TCC that is 
accurate to within a desired precision and can be evaluated rapidly using a finite number 
of arithmetic operations. 

In yet another aspect of the present invention, a computation of the partially 

20 coherent, defocused imagery of an arbitrary pattern in the paraxial approximation is 
provided in which the image intensity is represented purely in terms of analytical 
functions that is accurate to within a desired precision. 

The present invention further provides an expression for nonparaxial defocus 
imagery, in which an analytical expression for the nonparaxial defocus image is obtained 

25 by rescaling the coordinates of the paraxial defocused field. Again, this can be achieved 
to within a desired precision. Therefore, the inventive method can be extended to higher 
values of the NA with no additional computational effort and only a small, but known, 
increase in the overall error. 
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In another aspect of the present invention, aberrations are included by providing 
an analytical expression for an aberration pupil function by computing a Taylor 
expansion of the phase difference between the nonparaxial and paraxial phases. Each 
term in this expansion corresponds to a higher-order correction that is represented as a 
5 derivative of the paraxial result. With some additional computation depending on the 
precision desired, a pupil function including aberrations can be obtained by multiplying 
the inventive analytical paraxial defocus pupil function by the analytical aberration pupil 
function obtained using the Taylor expansion. This method is extendable to other phase 
errors that are described by the aberrations that may be present in a stepper. These 

10 aberrations are typically high-order polynomials in the frequency coordinate, and the 
coefficients are typically small - on the order of a few hundredths of a wavelength. 
Therefore, Taylor expansion is a natural way to extend this analytic methodology in the 
case of small aberrations. 

In yet another aspect of the present invention, the effects of amplitude variation 

1 5 across the pupil can be included. Such amplitude variations may become significant at 

higher values of the NA. These amplitude variations, or apodizations, come from energy 
conservation and vector effects that are not significant in the paraxial approximation. In 
accordance with the present invention, these apodization effects are also represented by 
approximated successively more accurately via an expansions, and the resulting 

20 amplitude effects are provided simply as multiplicative factors to the inventive paraxial 
defocused pupil function. 

The foregoing has outlined rather broadly the features and technical advantages of 
the present invention in order that the detailed description of the invention that follows 
may be better understood. Additional features and advantages of the invention will be 

25 described hereinafter which form the subject of the claims of the invention. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

For a more complete understanding of the present invention, and the advantages 
thereof, reference is now made to the following descriptions taken in conjunction with the 
accompanying drawings, in which: 
5 FIG. 1 is a schematic diagram illustrating a projection illumination (e.g. Kohler) system. 
FIG. 2 illustrates a diagram of source and pupil functions plotted in spatial frequency 
coordinates, including the TCC integration region. 

FIG. 3 is an illustration of a region of a mask pattern for which an image is to be 
simulated. 

1 0 FIGs. 4A-4D illustrate the four geometrical configurations for TCC integration regions 
for the case of a one-dimensional (ID) mask pattern. 

FIGs. 5A-5B illustrate a rotation of the coordinate system along the symmetry axis of the 
TCC integration region. 

FIGs. 6A-6R illustrate the 18 geometrical configurations for TCC integration regions for 
1 5 the case of a two-dimensional (2D) mask pattern. 

FIG. 7 illustrates the angles corresponding to the endpoints of the arcs defining the TCC 
integration region used in the contour integration in accordance with the present 
invention. 

FIG. 8 illustrates a plot of the error for a curve fit to the function -y/l-w 2 . 
20 FIG. 9 illustrates a plot of the error for an approximation of an integral used in computing 
TCCs according to the present invention. 

FIG. 10 illustrates a schematic diagram of a computer system capable of implementing 
the method of the present invention. 

FIG. 1 1 is a flow diagram illustrating a conventional method of computing an aerial 
25 image according to the Hopkins model. 

FIG. 12 illustrates a flow chart of a conventional method of computing TCCs. 

FIG. 13 illustrates a flow chart of an embodiment for computing TCCs in accordance 

with the present invention. 
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DETAILED DESCRIPTION 

In the following description, numerous specific details may be set forth to provide 
a thorough understanding of the present invention. However, it will be obvious to those 
skilled in the art that the present invention may be practiced without such specific details. 
5 In other instances, well-known features may have been shown in block diagram form in 
order not to obscure the present invention in unnecessary detail 

Refer now to the drawings wherein depicted elements are not necessarily shown 
to scale and wherein like or similar elements are designated by the same reference 
numeral through the several views. 
10 In accordance with the present invention, a fast, accurate analytical method for 

computing estimated aerial images, computed to within a desired tolerance. In 
accordance with the present invention, error estimates for the image computations are 
provided. The method of the present invention can incorporate aberrations. 

The source is assumed to be circular, preferably having a step function intensity. 
1 5 Alternatively, more general intensity distributions may be used, such as a Gaussian 

distribution. In addition, the pupil is assumed to be circular. The lens aberrations are 
assumed to be small (e.g. tens of milli waves), or zero. 

The aerial image intensity is determined according to the Hopkins model for an 
image projected by a given mask that is illuminated in a Kohler illumination (projection) 
20 system using a given set of illumination conditions. The illumination conditions are 

characterized by a source shape (S), a pupil function (P) and the mask pattern (M) 
expressed in spatial frequency coordinates. The Hopkins model assumes linear, shift- 
invariance of the optical system. This means two things: 1) (linearity) each distinct 
temporal frequency is unchanged by the system; and 2) (shift-invariance) each off-axis 
25 plane wave input into the system does not change the shape of the spectrum, but merely 
shifts it by some amount that depends on the angle of propagation. 

The following discussion is presented in order to provide a thorough 
understanding of the present invention, with specific examples presented as preferred 
embodiments. 
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1. General Theory 

Consider a modern microlithographic imaging system 200 (generally known as a 
Kohler illumination system) as illustrated in FIG. 1, comprised of a (finite) source 201, 
condenser optics 204, a semitransparent mask 207 to be imaged, objective optics 211 
having an entrance pupil 210, an exit pupil 212 and projection lens optics 218 (which 
typically includes multiple lenses) and a target, such as a silicon wafer in the image plane 
215. In an ideal system, the source 201 may be represented by point sources, each of 
which emits light rays, including on-axis rays 202 from the center of the source and off- 
axis rays emitted from an source point off of the optical axis 216. The on-axis rays 
correspond to a perfectly coherent source. A condenser lens 204 collects energy from the 
source and redirects the light rays (e.g. from the center of the source 205 and off-axis rays 
206) towards a mask 207 to be imaged. The on-axis and off-axis rays 208, 209 diffracted 
by the mask 207 are then collected by the entrance pupil 210 (i.e. the effective input 
plane) of the objective lens 21 1, is projected through the projection optics 218 and 
projects the mask image through the exit pupil 212 (i.e. the effective output plane) having 
an aperture 2 1 7 to form an aerial image at the target plane 215. Note that the aperture 
217 of the exit pupil 212, is typically a circle that filters out those light rays whose angle 
with respect to the optical axis 216 is larger than the quantity sin" 1 NA . The idea here is 
that each point of the source 201 is then mapped onto a unique plane wave, incident at 
some angle upon the mask 207. Further, a uniform source is mapped onto a uniform 
distribution of plane wave amplitudes upon the mask. The light that is transmitted 
through the projection optics 21 1 is partially coherent. 

The coherent electric field E at the point x (wherein the bold notation x and the 
notation x both refer equivalently to a vector quantity, hereinafter) on the wafer or image 
plane 215, due to the on-axis source point only, is given by the following expression: 




(1) 
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where M represents the spatial Fourier transform of the mask transmittance function, or 
frequency response of the mask 207 (i.e., the mask function), P represents the pupil 
function, or frequency response of the objective optics 211, NA is the numerical aperture 
of the objective lens at the exit pupil 217, and X is the illumination wavelength. 9? 2 
represents all of two-dimensional space. The vector quantity u represents a scaled 
frequency coordinate in which the pupil function P is zero when |u| > 1 . For an off-axis 

incident plane wave, the expression for the electric field due to a single source at the off- 
axis scaled source coordinate o becomes the following: 



The vector quantity NA <s represents the projection of a point on the unit sphere onto the 
mask plane 215, and in a Kohler illumination system, NA corresponds to the numerical 
aperture, and a represents a point in the illuminator 201, where the maximum magnitude 
of the vector o is equal to the ratio of the NA of the condenser lens to the NA of the 
projector lens. This ratio is called the partial coherence factor a 0 of the projection 
system 200. 

To simulate the effects of a partially coherent Kohler illumination system on the 
electric field intensity, the off-axis intensity is averaged over all of the plane waves 
generated by the illuminator. (Note that the use of partially coherent illumination 
represents both reality, as no illumination system can be perfectly spatially coherent, and 
desireability, as the shape of the illumination can be varied to achieve, say, optimum 
process latitude.) If it is assumed that the intensity is proportional to the square 
magnitude of the electric field distribution, then the intensity at a point x on the image 
plane 215 is modeled to take the form 



£(x,<t)= I d 2 uM(u-<y)P(u)exp -iln 



U X 



(2) 



jVcs(<^M 2 





(3) 
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where S is the source distribution for the illumination system. That is, as mentioned 
above, the intensity at the image plane is an average of the intensities emanating from the 
fields generated by the various point sources that comprise the illuminator. Substituting 
Eq. (2) into Eq. (3), the following expression for electric field intensity results: 



where the transmission cross-coefficient T (also referred to as the TCC or TCC integral) 
is given by 



In the integrand in Eq. (5), the two pupils functions P(u f + o) , ?(u" + o) will henceforth 
be referred to as the "offset pupils". The pupil function and source distribution function 



extent of the source and pupil apertures. This reduces the infinite integration interval 
specified in Eq. (5) into a finite interval, as shown in FIG. 2. 

FIG. 2 illustrates the regions of integration for the TCC integral of Equation (5), 
plotted in spatial frequency source coordinate space, where the coordinates (o x ,o y ), are 

relative to the center of the source function S(o) . Note that the source distribution 
function S(o) has a radius equal to the partial coherence factor <r 0 , and has a value that 
is zero outside of the boundary 404 where |o| > <r 0 . The two offset pupil functions 

P(u' + a) , P(u" + o) are represented by circles of unit radius 401 , 402 respectively, 
having centers offset from the center of the source function (referred to hereinafter, for 
convenience, interchangeably as the source) by the vector quantities -u f , -u", 
respectively. As discussed above, the offset pupil functions (referred to hereinafter, for 
convenience, interchageably as the pupils) take the value zero outside the boundaries 401, 
402. The two non-zero offset pupil functions (or pupils) have an intersection region 403. 
The region of integration 405 for the TCC integral of Equation (5) is the intersection of 





(5) 



has finite support; that is, |u| > 1 => P(u) = 0 and |u| > cr max => s(u)= 0, due to the finite 
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the nonzero source function 404 and the intersection 403 of the offset pupil functions 
401,402. 



10 



15 



20 



2. Periodic Objects 

As discussed above, in the Hopkins model, a simplifying assumption often made 
is that the mask 207, or a portion thereof, such as a local set of shapes, referred to as a 
cell, is considered to be a periodic object. The preferred embodiment of the present 
invention also makes use of this periodic assumption. Referring to FIG. 3, consider a 
semitransparent, periodic object 300 comprising cells 302, 303, 304 that are repetitions of 
a mask cell 301 with x- and y-periods p x and p y , respectively. Such an periodic mask 

300 may be represented mathematically as a Fourier series: 



M(u)= £ E C *» S 



( 

m 


X ' 




( 

n 


x) 


u x 




} 


Uy 

I Py 




V Px 


NA 







(6) 



where the c mn are the mask Fourier coefficients. The indices m and n will represent the 
various orders, or directions, into which the mask pattern 300 will split incoming light. 
Of course, the entire spectrum is not used in determining the image, as the objective lens 
21 1 filters out the higher orders. When the mask function M of Eq. (6) is substituted 
into the image intensity of Eq. (4), a series for the intensity results: 



Px Py 



7 ( X ) = X Z Z X ^viv^mVmVexP 
m'=^x> /w"=-oo //'--oo n"=-oo 

where the TCC integrals are evaluated at discrete lattice points due to the periodic 
assumption, and can be written as: 



,(7) 



1 ™' M '™' , M ' r — 1 

mnmn 



(8) 



{9 X <ly 1x Qy ) 

where the quantity q = pNA I X . The discrete lattice points, where the TCC integrals 
T mVmV are evaluated, are spaced at intervals of Au x = \jq x and Au v = \jq y . In 
principle, one could sum Eq. (7) directly, as the quantity T mWn , is bound to vanish when 
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the indices are large enough, making the sum finite. To see how this works, we take a 
closer look at the definition of the TCC of Eq. (5). 

The summation limits in Eq. (7) are now deduced from the bounds of the TCC 
integrals. Two geometrical conditions must be met simultaneously for nonzero values of 
the TCCs: 1) the separation between the centers of the two offset pupil functions must be 
less than twice their radius; and 2) both pupil functions must intersect the source 
function, which can be seen by reference again to FIG. 2. Letting k = m'-m an d 
l = n'-n\ these conditions (1) and (2) translate into Eq. (9) and (10), respectively, as 
follows: 



which provide finite integration limits. 

Note that, with this variable change, the intensity now takes the form of a Fourier 
series, which is expected from the periodic assumption. The conditions in Eqs. (9) and 
(10), and the change of variable suggested above, lead directly to the following exact 
finite sum for the intensity: 



u'-u1<2=> 




max 



ju'|, |u'|} < 1 + o => maxj + 




(10) 



M x M y {k) 





k=-M X l=-M y {k) 



< Px Py)\ 



(11) 



where 



M x =\2q x \ 



(12) 




(13) 



where the function \ x J is a floor function, representing the largest integer less than x (a real 



number), and 
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N x -k6{k) N y (m,k)-10(l) 

^kl = ^ ^ ^ ^ j c m+k,n+l c mn'^m+k,n+l,m,n > 

m=-ke{-k)-N x n=-W{-l)-Ny{m,k) 



9x 



(14) 

(15) 
(16) 



and 0is the Heaviside step function: 



1, x>0 
0, x < 0 ' 



(17) 



The quantities c m+k ll+l and c mn are the discrete coefficients corresponding to the mask 

functions M,M* in the expression for intensity l(x) of Eq. (4). The number of TCCs 
needed can be deduced from Eqs. (1 1) and (14); however, it is difficult to get an intuitive 
feel for how this number behaves as a function of the optical parameters. Eqs. (9) and 
(10) describe integer lattices bounded by ellipses, and therefore the number of points 
bounded by the ellipses is roughly equal to the product of the ellipses areas. An 
approximation to the number of TCCs, T, in the limit of a large unit cell area, is 



2NA Y 



PxPy 



(l + a)NA 



PxPy 



Note that this number is roughly proportional to the square of the period cell area. 

3. Characterization of the Integration Region 

Referring again to FIG. 2, the computation of the in-focus TCC is reduced to 
simply computing the common area of three circles, two of which are offset from the 
origin and have unit radius (the pupils), while the other is centered at the origin and has 
radius of a 0 (the source), Kintner (E.C. Kintner, "Method for the calculation of partially 

coherent imagery", Applied Optics 17, pp. 21 41 -21 S3 (1978)) showed that, for the case of 
one-dimensional (ID) masks (i.e. diffraction gratings) where the centers of the pupils are 
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constrained to lie on the horizontal (or vertical) axis, and where cKl, the TCC integration 
region geometry can be categorized into one of only four distinct configurations of the 
intersection 502 of the offset pupils 401, 402 and of the source 404, as illustrated in FIGs. 
4A-4D. These configurations are as follows: 
5 1) the source 404 is completely contained within the pupil intersection 502 to 

form an integration region 503 as in FIG. 4A; 

2) the source 404 protrudes out of one side of the pupil intersection 502 to form 
integration region 506 as in FIG. 4B; 

3) the source 404 protrudes out of both sides of the pupil intersection 502 to form 
10 integration region 509 as in FIG. 4C; and 

4) the pupil intersection 502 is completely contained within the source 404 to 
form an integration region 512 as in FIG. 4D. 

Clearly, the situation for the general two-dimensional (2D) mask is drastically 
more complicated. The number of distinct geometrical configurations for a 2D mask has 

1 5 not been published prior to the present application; the typical approach to computing the 

2D TCC has been to employ the minimum grid that contains the pupils and the source. 
Enumerating these integration regions is an important aid in determining the proper 
integration limits for analytical and adaptive 2D numerical integrations. It will be shown 
below, however, that, taking into account all possible symmetries, there are only 1 8 

20 distinct geometrical configurations to consider. 

In considering the 2D case, note that the coordinate system can always be rotated 
such that the symmetry axis of the intersection region of the offset pupils is vertical (that 
is, aligned along one of the source coordinate axes). For example, FIG. 5A illustrates a 
source coordinate system having horizontal axis a x and vertical axis a v centered on a 

25 source 404. Also illustrated is an offset pupil function 401 having center coordinates 

(w',v') offset from the center of the source 404 and a second pupil 402 offset at (u\v n ) 
relative to the source 404. The pupils 401, 402 have an intersection region 502 whose 
axis of symmetry 601 is tilted with respect to the vertical axis cr y . A coordinate rotation 
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can be performed so that the new rotated vertical axis o' y is parallel to the axis of 
symmetry 601 of the pupil intersection region 502, as illustrated in FIG. 5B. The pupils 
401, 402 will have new center coordinates (w',v), (w",v)> respectively. Note in particular 
that both pupil functions 401, 402 now have the same vertical center coordinate v in the 
5 rotated coordinate system [a' x , & y ) . By performing this rotation, the myriad of possible 

geometrical combinations is limited, and only a single vertical coordinate is necessary to 
specify the TCC. Organizing the different cases is also simplified, as the 2D TCC can 
now be viewed as a perturbation of the ID TCC. The various geometries can then be 
organized by their configuration with the y-components of the frequencies set to zero. 
10 Before enumerating the different cases, the desired change of coordinates is 

stated. Let the tangent of the rotation angle be denoted as 



25 



v -v 
r- 



u -id- 
Then the horizontal offsets and the common vertical offset are as follows: 

u' + rv' u+rv" 



Vl + r 2 5 VTf 



r 1 



(18) 



(19a) 



15 v^. (19b) 

Vl + r 2 



Also, denote u x = max{u\u n }, u 2 =min {u\u n }, and v = v . 

In the following descriptions of the integration regions of interest for the two- 
dimensional case, the organization is as follows. The two-dimensional geometric 
configurations may be categorized by analogy to the four one-dimensional cases. In each 
20 of these four categories of geometries, the pupils are moved together vertically, and each 

change in geometry is noted. Note that one geometric configuration can be derived from 
more than one logical condition. The cases are illustrated in FIGs. 6A-6R, in the 
coordinate axes have been rotated around the center of the source function 404 to align 
with the axis of symmetry of the pupil intersection region 502. 



3.1) The Source is completely contained within the Pupil Intersection 
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This is achieved when u x - 1 < -a 0 and u 2 + 1 ><r 0 . 

In this first category, the changes in geometry as the pupil intersection region 502 
is moved vertically relative to the source function 404 are considered. It is assumed that 
\u { | > |a/ 2 1 ; the complementary case results in the mirror image about the vertical The 
following quantities are defined: 

v, =V(l-^o) 2 -"i> v 2 = Vl-fo+"i) 2 . (20), (21), 

v 3 = ^(l-o- 0 ) 2 -u 2 2 , v 4 = ^\-(a 0 -u 2 ) 2 . (22), (23), 



= ^-^i-»iY . v, = ^ 0 2 -^("i+w 2 ) 2 , (24), (25), 

v 5 = v a ~ v /? . v 6 = v a + v /? . (26), (27), 

10 v 7 = ^(1 + cTo) 2 -^ 2 , v 8 =V( 1 + ^o) 2 -w 2 , (28), (29), 

V9=V^o 2 -(l-^) 2 , v IO =V CT o 2 -(l-" 2 ) 2 5 (30),(31), 
v„ = V^o 2 -(l + "i) 2 > v l2 =V< 7 o 2 -( 1 + ^) 2 > (32), (33), 

v i3 = ^-(^o-hl) 2 . v 14 = ^l-(cr 0 -|w 2 |) 2 . (34), (35). 

Then the distinct geometrical cases are determined under the following conditions: 
15 Case 1 (FIG. 6A): |v|<v,. 

Case 2 (FIG. 6B): w 2 < 0 , v 2 < v 3 and v, < |v| < v 2 . 

u 2 < 0 , v 3 < v 2 and v, < |v| < v 3 . 

u 2 > 0 , v 2 < v 5 and v, < |v| < v 2 . 

u 2 > 0 , v 2 > v 5 and v x < jv| < v 5 . 

20 Case 3 (FIG. 6C): * 2 < 0 , " 2 < v * and v 2 < |v| < v 3 . 

u 2 > 0 , v 2 < v 5 and v 2 < |vj < v 5 . 
Case 4 (FIG. 6D): u 2 < 0, v 2 < v 3 , v 4 < v 5 , and v 3 < |v| < v 4 . 
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u 2 < 0 , v 2 < v 3 , v 4 > v 5 , and v 3 < |v| < v 5 . 

u 2 < 0 , v 2 > v 3 , v 4 < v 5 , and v 2 < |v| < v 4 . 

u 2 < 0 , v 2 > v 3 , v 4 > v s , and v 2 < |v| < v 5 . 
Case 5 (FIG. 6E): u 2 < 0 , v 4 < v 5 , and v 4 < |v| < v 5 . 
Case 6 (FIG. 6F): u 2 < 0 , v 4 < v 5 , and v 5 < |v| < v 6 . 

v 4 >v 5 , and v 4 <|v|< v 6 . 

Case 7 (FIG. 6G): u 2 < -^^-m, , and v 6 < Ivl < v 7 . 

1 + (T 0 

Case 8 (FIG. 6H): u 2 < 0 , v 4 > v s , and v 5 < |v| < v 4 . 
Case 9 (FIG. 61): u 2 > 0 , v 2 > v 5 and v 5 < |v| < v 2 . 
Case 1 0 (FIG. 6J): u 2 < 0 , v 2 > v 3 , and v 3 < |v| < v 2 . 



3.2) The Source protrudes out of one side of the Pupil Intersection 

Recall that only the case of |w, | > \u 2 \ is being considered; the other case follows 
by symmetry. This is achieved when m, — 1 > -cr 0 and w 2 + 1 ><r 0 , and the distinct 
geometrical cases are determined under the following conditions: 

Case 1 1 (FIG. 6K): u 2 < - ] —^-u x , v 9 < v 3 , |v| < v 9 . 

1+(T 0 

u 2 < - - — u x , v 9 > v 3 ,|v| < v 3 . 
l + cr 0 

l -cr o ™ II 

— W| <u 2 < 2(7 0 — u, , v<v 9 . 

l+<7 0 



M 2 >2C7 0 -M,, |V|<V 9 . 

Case 12 (FIG. 6L): w 2 < -— ^-m, , v 9 > v 3 , v 9 < v 10 , v 5 < v 9 , v 3 < |v| < v 5 . 

1 + <T 0 
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m 2 <-^— — v 9 >v 3 , v 9 <v 10 , v 5 >v 9 , v 3 <|v|<v 9 . 
l + cr n 



u 2 <-\-^-u x , v 9 >v 3 , v 10 <v 9 , v 5 >v 10) v 3 <|v|<v 10 . 
l + cr 0 

u 2 <-\— — v 9 >v 3 , v 10 <v 9 , v 5 <v 10 , v 3 <|v|<v 5 . 
l + cr 0 

Case 13 (FIG. 6M): w 2 < -^^-Wi , v 9 > v 3 , v 10 < v 9 , v 5 > v I0 , v 5 > v 9 , v 3 < |v| < v 9 . 

l + cr 0 

u 2 < > v 9 >v 3 , v 10 <v 9 , v 5 >v 10 ,v 5 <v 9 , v 3 <|v|<v 5 . 

l + cr 0 

Case 14 (FIG. 6N): u 2 < -— ^-w, , v 9 > v 3 , v 10 < v 9 , v 5 > v 10 , v 5 < v 9 , v 5 < |v| < v 9 . 

l + cr 0 

u 2 <-j— v 9 >v 3 , v 10 <v 9 , v 5 >v l0 ,v 5 <v 9 , v 5 <|v|<v 9 . 

l+Cfr, 



l + <T ft 



W| <w 2 <2cr 0 -« 1 ,v 5 <v 9 ,v 9 >v 10 ,v 10 <|v|<v 9 . 



10 Case 15 (FIG. 60): u 2 < - j— — w, , v 9 > v 3 , v 9 < v 10 , v 5 < v 9 , v 5 < |v| < v 9 . 

l+<7 n 



l-O", 



u 2 < — °-u x , v 9 > v 3 , v 10 < v 9 , v 5 < v 10 , v 5 < |v| < v 10 . 



1+«T, 



1 + <T 0 

l-g 0 
l + (Tn 



«, <M 2 <2CT 0 -M,,V 5 <V 9 ,V 9 <V, 0 ,V 5 <|V|<V 9 



m, <w 2 <2tr 0 -M,,v 5 <v 9 ,v 9 >v 10 ,v 5 <|V|<V, 



"io- 



1 5 There are also many cases that revert back to a region defined in subsection 1 : 

l-cr n 



Case 3: 



u 2 < — -m, , v 9 < v 3 , v 9 < |v| < v 3 . 

l + <T n 
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1-cr, 



l+<7, 



<u 2 < 2cr 0 -u x , v 9 < v 5 , v 9 < |v| < v 5 . 



u 2 > 2cr 0 , v 9 < v| < v 13 . 



Case 4: 



Case 5: 



Case 6: 



u 2 <- 



u 2 < 



u 2 < 



u 2 <■ 



u 2 < 



u 2 < 



u 2 < 



u 2 <■ 



u 2 < 



u 2 < 



u 2 < 



u 2 <■ 



u 2 < 



~CT n 



+ <r n 



-cr n 



v 9 <v 3 ,v l0 <v 5 , v 3 <|v|<v, 



v 9 <v 3 ,v 10 >v 5 , v 3 < |v| < V, 



-cr n 



+ cr n 



-cr n 



-w, , v 9 > v 3 , v 9 < v l0 , v 5 > v 9 , v 10 < v 5 , V 



•"i , v 9 > v 3 , v 9 < v l0 , v 5 > v 9 , v 10 > V 5 , V 



+ cr n 



-u xt v 9 <v 3 ,v l0 <v 5 , v I0 <|v|<v 5 . 



+ <7 ft 



-w,, v 9 >v 3 , v 9 <v 10 , v 5 >v 9 , v 10 <v 5 ,v 



-cr n 



-<J« 



-"i , v 9 > v 3 , v 10 < v 9 , v 5 > v 10 , v 5 > v 9 , V 



-cr n 



"n V 9 <V 3 > V 10 <V S> V 5 <M <V 6- 



"l> V 9 <V 3 > V 5 <V 10> V 10 <| V | <V 6- 



-cr, 



+ <T, 



-cr n 



~f i , v 9 > v 3 , v 9 < v 10 , v 5 < v 9 , v 10 < |v| < 



-«,, v 9 >v 3 , v 9 <v 10 , v 5 >v 9 , v 10 <v s , 



-cr, 



+ CT ft 



— Wj , v 9 >v 3 , v 9 <v 10 , v 5 >v 9 , v 10 >v 5 , 



~<7« 



+ cr n 



-M,, V 9 >V 3 , V |0 <V 9 , V 5 >V 10 ,V 5 >V 9 , 
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u 2 <-]— — v 9 >v 3 , v 10 <v 9 , v 5 >v 10 ,v 5 <v 9 , v 9 <|v|<v 6 . 

l + (7 0 

-— — «, <" 2 <2<7 0 -M,, V 9 <V 5 , V, 0 <|v|<V 6 . 
1 + ^0 

-i— <m 2 <2<7 0 -«,,v s <v 9 ,v 9 <v 10 ,v 10 <|v|<v 6 . 
1 °fi-K, <« 2 <2<T 0 -«,,v 5 <v 9 ,v 9 >v 10 ,v 9 <|v|<v 6 . 



l + <7 0 



1-<T 0 



Case 7: u 2 < -w, , v 9 > v 3 , v 9 < v 10 , v 5 < v 9 , v 9 < |v| < v 10 



l + <7 0 



1-<T 0 



w 2 <-- Mj > v 9 >v 3 , v 9 <v I0 , v 5 >v 9 , v 10 >v 5 ,v 5 < v <v 10 . 

\ + a 0 

— u, < w 2 < 2cr 0 - w, , v 9 < v 5 , v 6 < |v| < v 7 . 
l + cr 0 

- — w, <u 2 <2a 0 -u x ,v 5 <v 9 ,v 9 <v I0 ,v 9 <|v|<v 10 . 

l + <7 0 

— m, <w 2 <2(7 0 -« 1 ,Vj <v 9 ,v 9 <v l0 ,v 6 <|v|<v 7 . 

U 2 > 2<T 0 -M, , V, 3 < |v| < V 7 . 

1-C7 n I I 

Case 8: u 2 < -w, , v 9 < v 3 , v 10 < v 5 , v l0 < v < v 5 . 

~\— — "i < " 2 < 2 °"0 , V, < V 5 , V 5 < |v| < V, 0 . 

l + cr 0 



3.3) The Source protrudes out of both sides of the Pupil Intersection 

This is achieved when w, - 1 > -a 0 , u 2 + 1 < <r 0 , and cr 0 2 < w,w 2 + 1 , and the 
distinct geometrical cases are determined under the following conditions: 
Case 1 6 (FIG. 6P): v 5 < v l4 , Ivl < v 5 . 
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V 5 > V 14»| V |< V 14- 

Case 1 7 (FIG. 6Q): v 5 < v 14 , v 5 < |v| < v l4 . 

Cases that revert back to regions defined in subsections 1 and 2: 

Case 5: v 5 > v 14 , v 5 > v 13 , v l3 < |v| < v 5 . 

5 Case 6: v 5 < v 14 , v u <|v|< v 6 . 

v 5 >v 14 ,v 5 <v 13 ,v 13 <|v|<v 6 . 
v s >v 14 ,v 5 >v 13 ,v 5 <|v|<v 6 . 
Case 13: v 5 >v 14 ,v 5 <v 13 ,v 14 <|v|<v 5 . 

V 5 >V 14> V 5 >V B> V 14 <H <V 13- 

10 Case 14: v 5 < v 14 , v 14 < |v| < v 13 . 

V 5 > V 14 , V 5 <V, 3 ,V 5 < |v| < V, 3 . 

3.4) The Pupil Intersection is completely contained within the Source 

This is achieved when w, - 1 > -<r 0 , u 2 + 1 < cr 0 , and <x 0 2 < w,w 2 + 1 , and the 
distinct geometrical cases are determined under the following conditions: 
15 Case 18 (FIG. 6R): |v|<v 5 . 

Cases that revert back to regions defined in previous subsections: 
Case 6: v 13 < |v| < v 6 . 

Case 14: v i 4 < | v | <v i3- 
Case 17: v 5 <| v | <v i4- 



20 



4. Simplification via Contour Integration 

The main result of the previous subsection was to identify all of the possible 
geometrical forms for the regions of integration for the TCCs, and the conditions under 
which each particular form would occur. 
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Consider the case of in- focus imaging, where the image plane is coincident with 
the focal plane of the projection lens. The pupil function P(u) takes the value 1 when 
|u| < 1 and zero when |u| > 1 . Then the in- focus TCCs are simply the areas of their 

respective regions. Given that each of the bounding contours of these regions is a 
circular arc, the determination of the areas exactly in terms of analytic functions is 
straightforward. There would be 18 separate analytical expressions to be evaluated 
according to the integration regions defined above, which is unwieldy to implement. 
However, it is possible to simply the 18 analytical expressions into a common form, as 
disclosed in Gordon (see Gordon, R.L., "Exact Computation of Scalar, 2D Aerial 
Imagery," Proc. SPIE, Vol. 4692, pp. 517-528 (2002)). 

This is accomplished by using Green's Theorem to convert the area integrals to 
contour integrals. According to Green's Theorem, the expression of the area A (e.g., the 
TCC for the in- focus case) of a closed region D can be expressed in terms of its bounding 
contour: 

4Z))= jVn=~ j(<r x da y -cT y dcT x ). (36) 

D dD 

Because the closed contour consists of nothing more than a set of circular arcs, 
the parametrization of the above contour integral is achieved using the following 

parametrization: a x = u^ k) + p k cos<2>, a y = v^ k) + p k sin^ ; the final result takes the 

form 

where the following notation applies: 

p k is the radius, and {u^\ vffl) is the center of the circle containing the kth arc. 
<ffl and ffl are the angles subtended by the endpoints of the kth arc at its circle's center 
measured clockwise with respect to the positive horizontal axis a x . For example, FIG. 7 
illustrates a source function 404 and two offset pupil functions 401, 404 having center 
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coordinates 801, 802 respectively. The intersection of the two pupil functions 502 is the 
region enclosed by the two circular arcs having endpoints A, B. The integration region 
805 is defined by the intersection of the source function 404 and the two pupil functions 
401, 402, and is the region enclosed by the arcs having endpoints C, D and B. In this 
example, there are three circular arcs: 1) the arc having endpoints C, D having radius p x 
of the source function; 2) the arc having endpoints C, B having radius p 2 corresponding 
to the pupil function 402 centered at coordinate 802; and 3) the arc defined by endpoints 
B, D having radius p 3 corresponding to the first pupil function 401 centered at 801. The 

angles corresponding to the endpoints of arc C-D are respectively; for arc C-B 

the angles are $ 2 \ffl , respectively; and for arc C-D the angles are $ 3 \ffi , 
respectively, where all of the angles are measured clockwise with respect to the positive 
horizontal axis. In addition, sf^ = sin <f>^ , cf^ = cos <f>f\ j e {l,2} , and N is the number of 

arcs in the contour. 

There are two immediate advantages to the contour integration method: 

1) The number of unique geometrical configurations of the integration regions 
shrinks from 1 8 to 9, which would simplify the used used to compute TCCs. This 
number is less because some of the integration regions detailed above have the same type 
ofboundary. 

2) The final expression for the in- focus TCC takes the form of Eq. (37) which is 
evaluated with a single arctangent and a few square roots per arc, resulting in a more 
efficient algorithm. 

5. Extension to Defocus 

For scalar, moderate-to-high NA systems (i.e., NA between about 0.5 to 0.7), the 
pupil function of a defocused system (i.e. at a defocus position z along the optical axis) 
generally takes the form 




(38) 
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Generally speaking, one must attack the problem of computing the TCCs with the pupil 
function in Eq. (38) numerically. However, as previously discussed, numerical 
computations of TCCs are slow and have errors that are not well understood. In 
accordance with this invention, the TCCs for defocused systems (i.e. to compute 
defocused images) are expressed in terms of standard, analytical functions, as described 
in more detail below. This is to be done without sacrificing accuracy, as other 
algorithms, such as SOCS or bilinear decomposition, must do. 

A standard approximation for the defocused pupil function of Eq. (38) is the 
quadratic-phase, or paraxial, approximation, which is useful for small values of NA (e.g. 
less than about 0.5), which leads to the following expression for the defocused, paraxial 
pupil function: 

P^(u) = exp(^-/>^A^ 2 u 2 j. (39) 

Substituting P DP of Eq. (39) into the TCC integral of Eq. (5) leads to the following 
expression for the TCC: 

r(u',u") = C J^ 2 (rexp(-/2^a c), (40) 
DnS 

where D is the intersection of the pupils, S is the source, a = -NA 2 z(u f -iT)/ X , and C is a 



phase term, where C = exp -inNA 2 —{u' 2 -u nl ) . 

\ A J 



The TCC integral in Eq. (40) now has the problem of having a nonuniform 
integrand over the intersection region Dr\S . By analogy to the in- focus case discussed 
above, Stokes 1 Theorem can be applied to convert the TCC area integral of Eq. 40 into a 
contour integral. This is shown as follows. Consider the vector-valued function: 

VAF = exp(-/2^a-(r)z (41) 
where z is a unit vector in the direction of the normal to the o plane. Stokes Theorem 
converts the TCC double integral in Eq. (40) to a single, contour integral around the 
boundary of the integration region: 
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JV<t-VaF= jdo-F. (42) 
DnS d{DnS) 
The idea is to find the value of F , given V a F . F can be found exactly because it is a 
solution of the vector wave equation. That is, F can be written as 

F = — L^yt^f eX p(- i2/r a • o), (43) 
/2;r a x + fl v 

where x is a vector in the a x direction and y is a vector in the cr y . 

The expression in Eq. (43) is then put into the right-hand side of Eq. (42). The 
reduction of this contour integral into a simple, single integral will depend on the 
boundary arc. The boundary of the integration region is made up of anywhere between 1 
and 4 circular arcs with distinct centers and radii. Actually, there are only two distinct 
types of arcs: 

2 2 2 

1) the boundary of the source function in a-space is described by <j x + o y = cr 0 ; and 

2) the boundary of the pupil function in a-space is described by 

((T x -U x ) 2 +(<T y -U y f =1 

If these arcs are parameterized by the angle subtended from the center of their 
respective circles, then it turns out that both types of boundaries lead to the following 
single integral expression for the TCC, in accordance with the present invention : 

T{u\n') = - eXP 0 ( " I f ) t Pj Wn Pj |a|, <* 0 + ^))-/(2^|a|^ 0 + ^)J (44) 



i2;r|a| 
where 

[<x 0 , source 
Pj \U Pupil 



(45) 



= Tan" 1 ^, (46) 



&i*\#2 are ^ e ang^s corresponding to the endpoints of the arc, as illustrated in FIG. 7, 
and N is the number of circular arcs on the contour, and: 
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j(a,P) = \dfi'oM ft exp(- ia cos /?') (47) 
o 

Thus, the TCC in Eq. (44) is now expressed in terms of single integrals, rather 
than the double integrals of prior art TCC expressions. This integral is therefore more 
easily computed by numerical methods. However, it would be more desirable to compute 
these single integrals analytically to improve speed and accuracy. Unfortunately, the 
integral in Eq. (47) cannot be reduced to an exact analytical form, except for a few 
isolated special cases. More preferably, the integrand of Eq. (47) is replaced by a 
different function that differs from the original integrand by at most, some small number 
6 . This new integrand function will be chosen such that the integral can be expressed in 
terms of purely analytical functions. 

First, note that Eq. (47) can be exactly rewritten in the following form: 

cos/7 

j(a, ft) = sin ft exp(-iacosfi)+ia J dw-yjl-w 2 exp(- / aw) . (48) 

i 

Again, the integral on the right in Eq. (48) is not expressible in terms of analytical 
functions. Note in particular that + w . The first term vl - w has 

the property that it has a vertical tangent at w=l, but the integral of Vl - w has a known 



analytical solution. The second term V 1 + w is smoothly varying, but the integral does 
not have a known analytical solution. Therefore, in accordance with the present 



invention, an approximation for the integrand VI - w 2 is made that preserves the 



properties of the first term Vl- w (i.e. has a vertical tangent at w=l) and has smoothly 
varying properties, that is, the error decreases substantially monotonically as the number 
of terms in the approximation increases. In a preferred embodiment, the following 



polynomial approximation can be made for the integrand Vl - w 2 over the interval for w 
of [0,l] (which covers all possible values of cos /? ): 



Vl-w 2 =71^ 1 + £/>"' +7J. (49) 

V m=l J 
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For a given finite value of L, the coefficients f m may be determined by performing a 

curve fit of the polynomial to Vl-w 2 at L+l discrete values of w. In a preferred 
embodiment, for a given number of terms L, the coefficients f m are found by equating 
values of the left and right sides of Eq. (49) at equal steps over the interval for w of [0,l] 
and are deduced by solving a linear system of equations. It turns out that when L=4, the 
error \rj\ < 1(T 5 , which is sufficient for a typical single precision calculation. For 
example, referring to FIG. 8, a plot of the absolute value of rj 1101, i.e. the difference 

between the function Vl-w 2 on the left hand side of Eq. (49) and the right hand side of 
Eq. (49) for L=4 (i.e. the curve fit at points 1 106, 1 107, 1 108, 1 109 and 1 1 10). Note that 
the maximum error is about 2 x 1 (T 5 at the point w=0.08 ( 1 1 02). 

Substituting the polynomial expression of Eq. (49) into the integral of Eq. (48) 
leads to a series of integrals. To evaluate the integral of the first term of the polynomial 
approximation from Eq. (49) we can make use of the following integral expression: 

jdwy/]^cxp(-iaw^ ^a{\-w) +/exp(-/aw)a~ 1 Vl^w, (50) 

where 

F^) = \duJi^-A (51) 

o \ J 

is the Fresnel integral. Although the Fresnel integral is not typically recognized as a 
standard analytical function, it does, in fact, contain all of the properties of such 
functions, and, most importantly here, can be computed to within machine precision 
using a (countably) small number of arithmetic operations, for example, by using a 
continued fraction expansion that is highly convergent, as is well known in the art. 

The integrals from Eq. (48) including the higher order terms in the polynomial 
from Eq. (49) are computed using the relation 

~a~da) j dw ex P(~ iaw )- (52) 



FIS920030254US1 



34 



Eqs. (49), (50), and (52) are combined to realize, for L=4, the following approximation 
for the integral in Eq. (48) which is valid for w <= [0,l]: 

cos/? 

J JwVl-vw 2 exp(- iaw) = 
i 

a' 5 jo. 166667exp(- iau)\ll - w • 

[/2.o(-0.030585w 4 + 0.133654w 3 -0.358614w 2 + 1.49819w+3.o)a 4 + 
(-0.275266w 3 + 0.904993w 2 -1.69w+ 4.2390 l)ar 3 + (53) 
/(-0.963433w 2 + 2.12485w-2.22014)a 2 - 
(2.40858w-2.70556)a-/3.61287] 

-/0.443113exp(-/a)ar- 1/2 F[o.797585Va(l-w)}4.0a 4 -il.49847a 3 
+ 0.499087a 2 + / 0. 1 39995a - 1 .703 12)}+ e, 

5 where \s\ < 10" 5 for all values of a and we [0,l]. FIG. 9 illustrates a Log-Log plot of 
e vs a for various values of we[0,l]. The error curves 1102, 1103, 1104 and 1105 
correspond to the values of w where the maximum error r/ occurs in FIG. 8 at points 
1 102, 1 103, 1 104 and 1 105 on the w axis, respectively. Note that the various floating- 
point numbers in Eq. (53) come from the results of the curve fitting to obtain the 

10 coefficients f m in Eq. (49). Eq. (53) can be extended to values of w<0 by computing its 
complex conjugate. 

Therefore, the preferred embodiment in accordance with the present invention, the 
paraxially defocused TCC (i.e. for NA <0.5) of Eq. (44), i.e. : 

15 can now be computed analytically, within an error of 10" 5 , where the integrals J are 
expressed by the following expression, for L=4, as: 
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j(a,fi)= sin exp(- ia cos /?) 

+ iaa ~ 5 {o. 166667 exp(- iau}yj\- w- 

[/2.o(- 0.030585 w 4 + 0. 133654 w 3 - 0.358614 w 2 + 1.49819w + 3.o)fr 4 
+ (- 0.275266 w 3 + 0.904993 w 2 - 1 .69 w + 4.23901 V 

(54b) 

+ /(- 0.963433 w 2 + 2. 12485 w - 2.22014 )a 2 

- (2.40858 w - 2.70556 )cr - /3.61287 ] 

- i 0.443 1 13 exp(- /a >T 1/2 f[o.797585 Ja(l - w)](4.0a 4 - / 1 .49847 a 3 
+ 0.499087 a 2 + i 0. 1 39995 a - 1 .703 1 2 )} 

where a = 2;zp y |a|, a = -A^4 2 z(u'-u' r )//l, and >9 = ^ 0 +^ 2 (y) and <t> 0 +^ at the endpoints 

of the N arcs of the region of integration, which is the intersection of the source and 
offset pupil functions. Similar analytical expressions for TCC can be derived having 
greater accuracy by including additional terms in a substantially monotonic 
approximation of Eq. (48) (e.g. by using larger values of L in Eq. (49) or similar 
expansions of Eq. (48)). Thus, in accordance with the present invention, these TCCs can 
be obtained with a relatively small number of arithmetic operations relative to numerical 
integrations that are used in current simulators. Accordingly, the image intensities can 
likewise be computed with improved speed and accuracy. 
6. Extension to High NA 

Eq. (38) dictates the nonparaxial (i.e. where NA is greater than about 0.5) phase 
behavior of the defocused pupil function P D . Because of the square root in the 
exponential, a simple, analytical expression, such as that found in Eq. (53), which is the 
result of a paraxial approximation, is not possible. Nevertheless, there is a simple 
algebraic transformation of the electric field expression in Eq. (1), that depends only on 
the spatial (x) coordinate and not the spatial frequency, that corrects for the phase errors 
introduced by the paraxial approximation (see, for example, Forbes et al., "Algebraic 
corrections for paraxial wave fields," J. Opt. Soc. Am. A, Vol. 14, No. 12, pp. 3300-3315, 
Dec. 1997). By treating this nonparaxial correction as an algebraic transformation that is 
independent of spatial frequency, the values of the TCCs will not change. 
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The nonparaxial correction to the paraxial electric field intensity is given by the 
following transformation: Let I P (x,z) be the paraxial intensity, at a point x on the 
image plane, and z is the defocus amount. Then the nonparaxial intensity /(x,z) (i.e. 
for all NA including NA greater than about 0.7) can be approximated by 

/(x,z)=[l + A^ 2 ^^ (55a) 



x 



\/(NAz) 



g{x,z)= ■ (55b) 

1 + - 



<M 4 



256x*NA*z 4 

The factor 1 + NA 2 g(x, z) 2 in square brackets is analogous to an obliquity term that is 
absent in the paraxial approximation. The denominator on the right-hand side of Eq. 
(55b) corresponds to a switching function that turns off the correction in the focal plane. 
7. Extension to Aberrations 

A wave aberration may be represented as a phase term of polynomials, such as 
Zernike polynomials, across the pupil. For the purposes of lithography modeling in a 
production semiconductor manufacturing environment, these aberrations are typically 
small, for example, less than about 1-3% of the illumination wavelength X . It is 
convenient to express the aberrations as an exponential P A of a phase term across the 
pupil as: 



P,(u)=exp 



q=\ 



(56) 



where Z q is the q th Zernike polynomial, Q q is the corresponding Zernike coefficient, Q is 
typically 37 in most lithographic applications, and s w , which is a small, dimensionless 
number that represents the order of magnitude of the wavefront error (for example, as 
discussed above, e w may be on the order of 0.0 IX). Because these aberrations are 

considered small compared to the wavelength, the exponential pupil function P A can be 
Taylor expanded with respect to s w out to, say, to yield: 
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P,(u)«l+/2^ w £^(u)- 2*£^Z>) *£. (56a) 



9=1 L « =1 



Typically, aberrations are combined with a defocus term for simulation purposes 
in order to explore degradation of process window that a particular set of aberrations may 
cause. The combined pupil function is multiplicative; that is 3 it takes the form 
P(u) = P A (u)P D (u) . Aberrations in a lens will result in a lack of a single geometrical 
focus, and increases the size of the diffraction "blur" in the image. This degrades the 
ability to image critical features over a wide range of focus and dose conditions, which is 
what we mean by the degradation of the process window. In accordance with the present 
invention, two methods are presented for analytically computing the aberration pupil 
function P A . 

In a first embodiment, note that the phase term (i.e. n—NA 2 xi 2 ) of the paraxial defocus 

X 

pupil function of Eq. (39) is mathematically equivalent to the fourth Zernike polynomial Z 4 , 
more specifically, is a quadratic function of u, and thus can be considered as an effective 
aberration. That is, the paraxial defocus term can be combined with the <^ 4 aberration 
coefficient, thus containing the defocus amount z as a component. Substituting the series 
expansion of Eq. (56a) into the TCC expression in Eq. (5) results in a series of terms 
involving integrals of the form 



For example, if £2=8, then there are 21 such terms. Again, the double integrals defined by 
Eq. (57) can be simplified by using Stokes' Theorem (see Eq. (42)), into a single integral 
around the boundary of the TCC integration region: 



As in the previous case, the contour is made up of N circular arcs, the p arc having 




(57) 




(58) 




When the path is parametrized as follows: 
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10 



15 



20 



°x =U P +Pp C0S <f>> 

°x= V P +P P sin ^ 

with <f> e [0,2;r], and the angle <f> of the endpoint is measured relative to the horizontal 
axis, then the integral in Eq. (58) can be evaluated as a sum over functions defined by a 
single integral: 



(59) 



nm 2 



N m n 



p=\ k=0 £=0 







(I) 





uT k v n - e p k ; e x 



d(f> 



p 1 cos 2 <f> + p u cos<f> p 2 sin 2 <f> + p p v p sin^ 



(60) 



cos* ^ sin 



m + 1 h + 1 

It turns out that the integrals in Eq. (60) can be expressed in terms of incomplete beta 
functions: 

siir ^ I 



n 



, • K j l n f K + l L + l) _ 



a:+i i+i 



(61) 



where 



z 



(62) 



is the incomplete beta function. This function is well-known in the mathematical 
literature and algorithms exist that compute its value using only a few arithmetic 
operations. Such an expansion of the aberration pupil function in terms of incomplete 
beta functions has been considered previously for one-dimensional (ID) geometries (see, 
for example, Steel, W.H, "Effects of Small Aberrations on the Images of Partially 
Coherent Objects," J. Opt. Soc. Am., Vol. 47, No. 5, pp. 405-413 (1957)), but not for 
two-dimensional (2D) geometries. In accordance with the present invention, the form of 
Eq. (61) provides an analytical solution for the aberration pupil function P A for general 
2D geometries. This pupil function can be multiplied with the defocused pupil function 
P D , or for example P DP as provided in Eq. (39) of a preferred embodiment, and the 
resulting TCC computed, for example as in Eq. (54a), can be used to provide an 
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analytical expression for image intensity that includes aberrations, and that has small 
error. However, this solution is applicable to a limited range of defocus values, limited 
by the small value of s w that is assumed. 

More preferrably, in a second embodiment for computing aberrated images, the 
paraxial defocus term P DP is treated separately from the aberration term P A and the 
argument in the exponential of P DP is not assumed to be small (i.e. for z larger than the 
wavelength X ). In accordance with the present invention, the aberration pupil function 
P A is Taylor expanded to 2 nd order in s w in accordance with Eq. (56a), but the paraxial 
defocus pupil function P DP (see Eq. (39)) is kept as is. Substituting the Taylor expansion 
of P A into the TCC integral of Eq. (5), integrals of the form: 



The integral on the 2 nd line of Eq. (63) has the same form as the TCC integral of 
Eq. (40) and therefore, is expressed explicitly in terms of derivatives of the 
analytical functions of the form previously derived, for example as in Eqs. (54a) and 
(54b), whose integrals also have known analytical solutions. Therefore, the aberration 
pupil function P A is now expressed in terms of analytical functions which is then, in 
accordance with the present invention, multiplied with the paraxial defocused pupil 
function P DP to obtain an integrand for the TCC integral that includes aberrations. 
8. Energy conservation, vector effects, and air/resist interfaces 

In the simulation of the very high-NA (NA > 0.7) systems in use at the latest 
technology nodes, other physical effects such as obliquity due to power conservation 
constraints, vector diffraction effects, and the effects of the resist layer complicate the 
analytical structure of the Hopkins' TCC integral. These three effects, however, have the 
common property that they are essentially pupil apodizations; that is, they are amplitude 



j^ 2 a<<exp(-/2;ra.a) 



DnS 




(63) 



are obtained, where a = -NA 2 z(u' - u")/ X . 
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variations across the pupil. The functional forms that describe these variations each 
would make the TCC integral impossible to simplify in an exact form. For example, the 
propagation of the light into the resist from air produces the following pupil amplitude 
function: 

- p M 2V1-A14V 

5 M U ) = 1=T7 / , , , ' ( ( 

V1-M4V + J\-NA 2 vl 2 In 2 

where n is the index of refraction of the resist. Nevertheless, as with the aberration term, 
such apodizations can be considered as corrections to the unapodized image, so that the 
apodization terms can be Taylor expanded for small spatial frequencies. Thus, in a 
preferred embodiment of the present invention, an apodization pupil function P Ap , which 

10 can be expressed as an amplitude factor, would then be multiplied with the paraxial 
defocused pupil function P D , or P DP to obtain a pupil function that includes 
apodizations. 



9. Implementation in Computer Software 

15 The method of simulating the image intensity in accordance with the present 

invention may be implemented in machine readable code (i.e. software or computer 
program product) and performed on a computer, such as, but not limited to, the type 
illustrated in FIG. 10, wherein a central processing unit (CPU) 1300 is connected by 
means of a bus, network or other communications means 13 10 to storage devices such as 

20 disk 1301, or memory 1302, as well as connected to various input/output devices, such as 
a display monitor 1304, a mouse 1303, a keyboard 1305, and other I/O devices 1306 
capable of reading and writing data and instructions from machine readable media 1308 
such as tape or compact disk (CD). In addition, there may connections via a network 
13 10 to other computers and devices, such as may be exist a network of such devices, e.g. 

25 the Internet. Software to perform the method steps to obtain an image may be 

implemented as a program product and stored on machine readable media 1308. 
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FIG. 1 1 illustrates a typical flowchart of a method for computing images 
according to the Hopkins model. At step 101, characteristics of the imaging system are 
provided, such as the illumination wavelength X , the numerical aperture NA 5 the partial 
coherence factor a , the amount of defocus z and other parameters of the imaging 
5 system and process. A mask region to be imaged is provided (step 1 02) where typically 

the mask region is assumed to be periodic at a pitch p x9 p y in the x- and y-directions, 

respectively for a 2D mask. The spacing Au in the spatial frequency domain is 
determined by Aw = X/ pNA , which may be different in the x- and y-directions depending 
on the pitch p (see Eqs. (7) and (8) above). Next the summation limits from Eqs. (11) 

10 and (14) are determined (step 104). Other imaging parameters may be optionally 

provided, such as, for example lens aberrations (step 105). The TCCs of Eq. (5) are 
computed (step 106), typically by a numerical method to be described in more detail 
below. As an example, for a 40/um 2 mask cell, simulated at X = \ 92>nm , there are about 
3.6 million TCC integrals to compute. Mask parameters, for example in the form of 

15 polygon shapes, are provided (step 107). The Fourier spectrum of the image intensity for 

those polygons is computed (step 108), for example as disclosed in U.S. Patent 
Application Serial No. 10/353,900, filed on Jan. 28, 2003, which is commonly assigned 
to the Assignee of the present application. Optionally, the image can be simulated over a 
subset of the mask cell region, for example, one might be interested only in the CD of a 

20 line in a particular portion of the cell (step 109). Then the image intensity is computed 
over the simulation region according to Eqs. (14) and (11). 

More specifically, the conventional method for computing TCCs is illustrated in 
FIG. 12. First a pair of offset pupil functions are provided, by entering the center 
coordinates in scaled spatial frequency source coordinate space (step 901). Each offset 

25 pupil function represents a coupling between rays originating from different source points 

passing through the pupil. The region of integration (typically, between - 1 < o x < 1 and 

- 1 < a y < 1) is divided into an NxN grid block (step 902). For each grid block (step 

903), the location of the corners of each grid block with respect to the intersection regions 
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of the source function and the offset pupil functions are determined (step 905). Based on 
the location of the grid block corners, the area of the grid block within the integration 
region is determined (e.g. estimated) (step 906), and the TCC integral of Eq. (5) is 
computed by numerical summation of the source- weighted pupil functions over all grid 
5 blocks (907, 908), and scaled by the area of the integration region (step 909). Note that 
the accuracy of the integration by the prior art method is dependent on how fine the 
gridding is, but the computation time increases significantly as more grid blocks are used. 

A flow chart of an embodiment in accordance with the present invention is 
illustrated in FIG. 13. For each TCC, a pair of pupil functions having center coordinates 

10 offset in source coordinate space relative to the center of the source function (step 1001). 

The intersection region of the offset pupil and source functions are identified as the 
region of integration and classified according to one of the 18 geometries discussed in 
subsection 3 above (step 1002). The boundary arcs of the intersection region are then 
characterized in terms of their intersection points (i.e. the arc endpoints) and the angular 

1 5 coordinates corresponding to the arc endpoints are determined (step 1003). Next, the 

form of the TCC integrand is determined by approximating the pupil function by at least 
a paraxial defocus term times a polynomial in the a x ,cr y coordinates (step 1004). The 

contour integrals for the expansion of the pupil functions are determined (for example, 
using Eq. (54b)) (step 1005), and then the contributions from all arcs are added (step 

20 1006) and the TCC, in a form such as Eq. (54a), are computed (step 1007). If aberrations 

are to be included, an aberration pupil function is preferably obtained by an Taylor 
expansion of a deviation from a spherical lens, in the form of analytical functions of the 
form in Eq. (63) and used to modify the integrand in step 1004. Additionally, step 1004 
could be further modified to include apodization pupil functions as additional factors. 

25 Although the present invention and its advantages have been described in detail, it 

should be understood that various changes, substitutions and alterations can be made 
herein without departing from the spirit and scope of the invention as defined by the 
appended claims. 
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